cd /home/zhaoshuo/pan_niger
##ASM285v2 is the reference genome. Mapping all assemblies to this reference genome.

#path to the reference genome
ref_gnm="/home/zhaoshuo/pan_niger/ref_ASM1928827v1_genomic.fna"
#path to the assemblies
ASM285v2_gnm="/home/zhaoshuo/pan_niger/ASM285v2_genomic.fna"
ASM2576888v1_gnm="/home/zhaoshuo/pan_niger/ASM2576888v1_genomic.fna"
ASM2978390v1_gnm="/home/zhaoshuo/pan_niger/ASM2978390v1_genomic.fna"
ASM2978392v1_gnm="/home/zhaoshuo/pan_niger/ASM2978392v1_genomic.fna"
Asplac1_gnm="/home/zhaoshuo/pan_niger/Asplac1_genomic.fna"
CSR3_gnm="/home/zhaoshuo/pan_niger/CSR3_genomic.fna"
#mapping all assemblies to the reference genome
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o ASM285v2_to_ref.paf $ref_gnm $ASM285v2_gnm
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o ASM2576888v1_to_ref.paf $ref_gnm $ASM2576888v1_gnm
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o ASM2978390v1_to_ref.paf $ref_gnm $ASM2978390v1_gnm
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o ASM2978392v1_to_ref.paf $ref_gnm $ASM2978392v1_gnm
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o Asplac1_to_ref.paf $ref_gnm $Asplac1_gnm
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o CSR3_to_ref.paf $ref_gnm $CSR3_gnm

wc - l *.paf
#2769 ASM2576888v1_to_ref.paf
#1228 ASM285v2_to_ref.paf
#2417 ASM2978390v1_to_ref.paf
#1439 ASM2978392v1_to_ref.paf
#2082 Asplac1_to_ref.paf
#1883 CSR3_to_ref.paf
#11818 total

##generate a "genome file" for "bedtools complement" using infoseq
infoseq -length -sequence $ASM285v2_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > ASM285v2.gnminfo 
infoseq -length -sequence $ASM2576888v1_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > ASM2576888v1.gnminfo
infoseq -length -sequence $ASM2978390v1_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > ASM2978390v1.gnminfo
infoseq -length -sequence $ASM2978392v1_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > ASM2978392v1.gnminfo
infoseq -length -sequence $Asplac1_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > Asplac1.gnminfo
infoseq -length -sequence $CSR3_gnm | sed 's/>[^ ]*/>&/' | awk 'NR>1 {print $3"\t"$6}' > CSR3.gnminfo
#extract the alignments with length >= 300bp
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' ASM285v2_to_ref.paf | bedtools sort -g ASM285v2.gnminfo > ASM285v2_to_ref.align.bed
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' ASM2576888v1_to_ref.paf | bedtools sort -g ASM2576888v1.gnminfo > ASM2576888v1_to_ref.align.bed
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' ASM2978390v1_to_ref.paf | bedtools sort -g ASM2978390v1.gnminfo > ASM2978390v1_to_ref.align.bed
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' ASM2978392v1_to_ref.paf | bedtools sort -g ASM2978392v1.gnminfo > ASM2978392v1_to_ref.align.bed
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' Asplac1_to_ref.paf | bedtools sort -g Asplac1.gnminfo > Asplac1_to_ref.align.bed
awk '{len=$4-$3;if(len>=300){print $1"\t"$3"\t"$4}}' CSR3_to_ref.paf | bedtools sort -g CSR3.gnminfo > CSR3_to_ref.align.bed

wc -l *.align.bed
#1283 ASM2576888v1_to_ref.align.bed
#767 ASM285v2_to_ref.align.bed
#2005 ASM2978390v1_to_ref.align.bed
#959 ASM2978392v1_to_ref.align.bed
#966 Asplac1_to_ref.align.bed
#997 CSR3_to_ref.align.bed
#6977 total

#get unaligned regions in the reference genome(>=500bp)
bedtools complement -i ASM285v2_to_ref.align.bed -g ASM285v2.gnminfo | awk '{if($3-$2>=500){print}}' > ASM285v2_to_ref.unalign.bed
bedtools complement -i ASM2576888v1_to_ref.align.bed -g ASM2576888v1.gnminfo | awk '{if($3-$2>=500){print}}' > ASM2576888v1_to_ref.unalign.bed
bedtools complement -i ASM2978390v1_to_ref.align.bed -g ASM2978390v1.gnminfo | awk '{if($3-$2>=500){print}}' > ASM2978390v1_to_ref.unalign.bed
bedtools complement -i ASM2978392v1_to_ref.align.bed -g ASM2978392v1.gnminfo | awk '{if($3-$2>=500){print}}' > ASM2978392v1_to_ref.unalign.bed
bedtools complement -i Asplac1_to_ref.align.bed -g Asplac1.gnminfo | awk '{if($3-$2>=500){print}}' > Asplac1_to_ref.unalign.bed
bedtools complement -i CSR3_to_ref.align.bed -g CSR3.gnminfo | awk '{if($3-$2>=500){print}}' > CSR3_to_ref.unalign.bed

wc -l *.unalign.bed
#787 ASM2576888v1_to_ref.unalign.bed
#57 ASM285v2_to_ref.unalign.bed
#305 ASM2978390v1_to_ref.unalign.bed
#175 ASM2978392v1_to_ref.unalign.bed
#82 Asplac1_to_ref.unalign.bed
#82 CSR3_to_ref.unalign.bed
#1488 total

#get unaligned fasta sequences
bedtools getfasta -name+ -fo ASM285v2_to_ref.unalign.fa -fi $ASM285v2_gnm -bed ASM285v2_to_ref.unalign.bed
bedtools getfasta -name+ -fo ASM2576888v1_to_ref.unalign.fa -fi $ASM2576888v1_gnm -bed ASM2576888v1_to_ref.unalign.bed
bedtools getfasta -name+ -fo ASM2978390v1_to_ref.unalign.fa -fi $ASM2978390v1_gnm -bed ASM2978390v1_to_ref.unalign.bed
bedtools getfasta -name+ -fo ASM2978392v1_to_ref.unalign.fa -fi $ASM2978392v1_gnm -bed ASM2978392v1_to_ref.unalign.bed
bedtools getfasta -name+ -fo Asplac1_to_ref.unalign.fa -fi $Asplac1_gnm -bed Asplac1_to_ref.unalign.bed
bedtools getfasta -name+ -fo CSR3_to_ref.unalign.fa -fi $CSR3_gnm -bed CSR3_to_ref.unalign.bed
sed -i 's/^>::/>/' *.unalign.fa

ll *.unalign.fa | awk '{print $5,$9}' 
#4.2M ASM2576888v1_to_ref.unalign.fa
#345K ASM285v2_to_ref.unalign.fa
#2.1M ASM2978390v1_to_ref.unalign.fa
#1.1M ASM2978392v1_to_ref.unalign.fa
#495K Asplac1_to_ref.unalign.fa
#506K CSR3_to_ref.unalign.fa

#combine all unaligned sequenses and cluster them 
cat *.unalign.fa > all.unalign.fa
source /home/SOFTWARE/mmseqs/util/bash-completion.sh
/home/SOFTWARE/mmseqs/bin/mmseqs easy-cluster all.unalign.fa all.unalign.mmseqs.cluster tmp
sed -i -e '/^>/s/:/_/' -e '/^>/s/-/_/' all.unalign.mmseqs.cluster_rep_seq.fasta
grep -c '^>' *.unalign.fa all.unalign.mmseqs.cluster_rep_seq.fasta
#all.unalign.fa:1488
#ASM2576888v1_to_ref.unalign.fa:787
#ASM285v2_to_ref.unalign.fa:57
#ASM2978390v1_to_ref.unalign.fa:305
#ASM2978392v1_to_ref.unalign.fa:175
#Asplac1_to_ref.unalign.fa:82
#CSR3_to_ref.unalign.fa:82
#all.unalign.mmseqs.cluster_rep_seq.fasta:1256
ll all.unalign.mmseqs.cluster_rep_seq.fasta | awk '{print $5,$9}' 
#7.5M all.unalign.mmseqs.cluster_rep_seq.fasta

# align representing unaligned sequences to MR1 genome to remove redundance.
minimap2 -t 110 -x asm10 -p 0.01 -N 100 -o allunalignUniq_to_ref.paf $ref_gnm all.unalign.mmseqs.cluster_rep_seq.fasta

#merge overlapping alignments
cut -f 1-4 allunalignUniq_to_ref.paf | sed -r 's/\t/_/' | sort -k1,1 -k2n,2 > allunalignUniq_to_ref.bed
bedtools merge -d 300 -i allunalignUniq_to_ref.bed | sed -r 's/_/\t/3' > allunalignUniq_to_ref.merged.bed4

#output sequences with overlapping coverage < 0.8 with reference genome
awk '{if($1!=lastseq){if(lastseq!="" && matches/lastlen<0.8 && lastlen-matches>500){print lastseq"\t"lastlen"\t"matches};lastseq=$1;matches=$4-$3;lastlen=$2}else{matches=matches+$4-$3}}END{if(matches/lastlen<0.8 && lastlen-matches>500){print lastseq"\t"lastlen"\t"matches}}' allunalignUniq_to_ref.merged.bed4 > allunalignUniq_to_ref.merged.nool
awk '{len=len+$2; matches=matches+$3}END{print NR"\t"len"\t"matches}' allunalignUniq_to_ref.merged.nool
#seq    len    matches
#288	2591099	68441

#retrieve these sequences
cut -f 1 allunalignUniq_to_ref.merged.nool > allunalignUniq_to_ref.merged.nool.names
sed -r -i 's/\s+$//' all.unalign.mmseqs.cluster_rep_seq.fasta
awk -F'>' 'NR==FNR{ids[$0]=1; next} NF>1{f=ids[$2]} f' allunalignUniq_to_ref.merged.nool.names all.unalign.mmseqs.cluster_rep_seq.fasta > allunalignUniq_to_ref.merged.nool.fa

#combine  reference genome and unaligned sequences
cat $ref_gnm allunalignUniq_to_ref.merged.nool.fa > niger_pan.fa










